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ABSTRACT 



The waveguide mode tropospheric propagation effect prediction program, M- 
Layer, originally written by Naval Command Control and Ocean Surveillance Center, 
Research, Development, Test and Engineering Division (NRaD), is revised for 
greater accuracy, speed and stability. The accuracy improvement is achieved first by 
converting the extended complex number representation into the representation by 
the complex exponent then by re-writing the group of Airy function computation 
subroutines. This accuracy improvement makes it possible to implement a self- 
checking procedure for determining the proper method to evaluate the height gain 
function. Finally, a new mode locating algorithm is introduced which improves the 
efficiency of mode search and eliminates the looping problem observed. The revision 
has been documented and the new program source code has been delivered to 
NRaD. It is also recommended that the mode search protocol, not just the mode 
locating algorithm introduced in this revision, be completely revised for better 
performance. 
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I. INTRODUCTION 



M-Layer is a FORTRAN program for computing the propagation factor of an 
electromagnetic (EM) wave in a stratified atmosphere. It is desirable to extend the 
capability of this program to include a layer of random medium representing the air- 
ocean interface. To achieve this goal, there are many basic theoretical problems 
which have to be answered. First of all, the effect of the earth curvature in this 
program is taken care of through the classical earth-flattening approximation [Ref. 1], 
but the result [Ref. 2] does not agree with the more recent diffraction theory of Fock 
[Ref. 3] near the surface of the earth. Then there is the question about the better 
method to model the atmospheric refractive index profile, either piecewise linear or 
quadratic, to be resolved by a new earth-flattening approximation under development 
at NPS. The new approximation will also determine the functions to be used for the 
representation of the EM fields in each layer through uniform asymptotic theories. 
Within some proper region, these new functions are expected to reduce to the Airy 
functions utilized by M-Layer. The evolutionary nature of this effort prompted this 
review to improve the inner workings of the M-Layer program. In particular, the 
subroutines to search for the modes and those for evaluating the Airy functions will 
remain as an important part of a program investigating questions about EM wave 
propagation by solving the related boundary value problem. 
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It can never be overemphasized that a boundary value problem which includes 
a layer of random medium or some range dependent inhomogeneity, set up according 
to the Maxwell equations, will include backscattering in its solution. This is in sharp 
contrast to those numerical procedures based on the parabolic approximation to the 
wave equation for which the backscattering is completely ignored. 

In what follows, the M-Layer program and the reasons for replacing the 
extended complex numbers with their complex exponent representations are 
discussed, together with some other problems encountered and resolved during this 
investigation. 

A. M-LAYER 

In M-Layer, the index of refraction of the atmosphere is assumed to be height 
dependent and is approximated with a continuous piecewise linear profile. The 
classical earth-flattening approximation is utilized to allow the use of the cylindrical 
coordinate system while retaining the effect of the curvature of the earth. This is 
done simply by substituting the index of refraction with the modified index of 
refraction, which also has a piecewise linear profile [Ref. 1]. 

The source of the EM radiation is assumed to be either a vertical electric 
dipole or a vertical "magnetic dipole", with the latter providing an approximation to 
the radiation of a horizontal electric dipole. The dipole is located along the positive 
z-axis of the cylindrical coordinate system while the origin is sitting on the ground. 
The x-y plane is the ‘flattened’ earth surface. After carrying out the Hankel 
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transform along the radial direction, the resulting spectrum of the Hertzian dipole 
field within each layer of a linear segment of the modified refractive index profile is 
reduced to a linear combination of the Airy functions. Specifically, the layers are 
numbered to increase with height, with the first layer being the one above the 
ground. The spectrum of the Hertzian dipole field is proportional to the product of 
the values, at the transmitter height and at the receiver height respectively, of the 
height-gain function. At a height within the i-th layer, the height-gain function is 
given by [Ref. 4]: 

, ( 1 ) 

where p is the radial component of the propagation vector and is also the spectral 
variable of the Hankel transform; hence it is the same throughout all layers. It is a 
complex variable whose imaginary part represents the radial attenuation rate of the 
spectral component of the Hertzian dipole field. Under the classical earth-flattening 
approximation, the spectrum of the Hertzian dipole field contains a discrete portion 
and a branch cut. The discrete spectrum gives rise to the creeping wave modes 
diffracted by the earth surface and the dielectric waveguide modes supported by the 
layered atmosphere. The contribution from the branch cut is usually negligible, 
especially for the field in the shadow of the earth. The M-Layer program locates the 
discrete spectrum for modes having a radial attenuation rate below a predetermined 
value. Contributions from these modes determine the propagation factor of the 
wave. 
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The variable q i in the i-th layer is a dimensionless linear function of height z 
with the free space wavenumber k, the modified index of refraction m i at the lower 
boundary z = z- v the slope of the modified index of refraction a;/2 and p as 
parameters: 
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The height dependence of the field is given in terms of the functions k x (q^) and 
k 2 (qj), which are proportional to the Airy functions Ai(-q p ]2K//3 ) and Ai(-q { ) 
respectively. Of these two functions, at a height so large that q i is large and positive, 
k^qj) represents a downward going wave and e^ 4 n ^ 3 k ^ (q ,) +k 2 (q ;) represents an 
upward going wave. The coefficients ^ i and5,- are determined by the conditions on 
the continuity of the Hertzian dipole field and its derivative across layer boundaries 
and by the normalization condition that the integral of the square of the height-gain 
function over all height equals unity. 

To fulfill the radiation condition, the highest layer is given the same refractive 
index as the free space above it and only the outgoing wave is allowed within this 
layer. Below the ‘flattened’ earth surface, the field is assumed to be a plane wave 
propagating downward. Hence, only the normalization factors are required in the 
highest layer and in the ground. By assigning#; to unity in the highest layer, all the 
coefficients A; and#,- can be determined, according to the boundary conditions, to 
within a multiplicative factor for#;. This multiplicative factor is then deduced from 
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the normalization condition. This procedure can also be carried out from the ground 
level up. That these coefficients can be computed either from the highest level down 
or from the lowest level up is a result of the fact that p belongs to the discrete 
spectrum of the Hertzian dipole field. Consequently, agreement between these two 
ways of evaluating theA ; - and£ ; - coefficients confirms that a mode has been located 
accurately. 

B. EXTENDED COMPLEX NUMBER REPRESENTATION 

The discrete spectrum of the Hertzian dipole field corresponds to the zeroes 
of the modal function which is a determinant whose elements consist of k^qf) and 
kjiqt) at the layer boundaries. Numerically, the magnitude of this modal function 
causes overflow and underflow problems as k^q-) or becomes exponentially 

large or small for complex q t values. In the M-Layer program, to overcome this 
problem, a complex number is written as a scaled number, which is complex, 
multiplied by a scaling factor which is an integer power of e, the base of natural 
logarithm. This integer is chosen so that the greater of the absolute values of the 
real part and the imaginary part of the scaled number lies within e 11 . A complex 
number written in this form is called an extended complex number. Multiplication 
of two extended complex numbers requires summing the two integer exponents in 
addition to carrying out the regular complex multiplication of the scaled numbers. 
Addition of two such numbers is achieved through the use of an addition subroutine: 
the larger scaling factor is factored out of both addends before they are combined. 
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The scaling factor is adjusted after each addition and after a sequence of 
multiplications to make sure that the resulting scaled number is still within the 
desired range. Addition is troublesome when the two numbers to be added nearly 
cancel each other. Under this circumstance, the scaling factors of the two numbers 
are identical and both the real parts and the imaginary parts of the scaled numbers 
are almost equal with opposite signs. It is clear that the real part and the imaginary 
part of the sum lose their accuracies to different degrees; hence the phase angle may 
incur substantial error. To remedy this situation, interpolation procedures have to 
be devised. 

As two complex numbers come close to cancel each other, they must be out of 
phase by almost 180 degrees. By factoring out the square root of their product 
instead of the scale factor, the resulting addends become reciprocal to each other, 
both lying within an identical small angle to, and on the same side of, the imaginary 
axis. They are close to the unit circle, but one is on the inside and the other is on 
the outside. Taking out further a phase factor of 7t/2 after writing the addends in 
their exponential forms, the exponents become small numbers for which a Taylor 
series expansion of the exponential function converges rapidly and can be used for 
interpolating the sum to achieve higher accuracy. Note that after the extra phase 
factor of 7t/2 is removed from the addends, it is actually the difference of the 
resulting two reciprocals which is computed. This procedure effectively picks the 
direction on the complex plane along which the addends are almost opposing each 
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other to carry out their cancellation. The resulting sum has a phase angle nearly 
perpendicular to this chosen direction. 

It is evident that the representation of a complex number by its complex 
exponent of base e provides better phase accuracy for addition. A one-to-one 
correspondence can be achieved by restricting the imaginary part of this exponent to 
within -n and t r. This will be called the exponential representation or the complex 
exponent representation henceforth. It is convenient for multiplication: adding the 
complex exponents of the two factors will suffice. Conversion of the M-Layer 
program from the extended complex number to the complex exponent representation 
has been carried out. 



C. OTHER REVISIONS 

As better precision is achieved, problems with the mode search procedure and 
the evaluation of the and B i coefficients become severe. They are thoroughly 
investigated and resolved. For mode search, although the division of the region of 
interest into "contour rectangles" and further into square "meshes", and the search 
pattern to r the sides of a "contour rectangle" to find and follow "phase 

lines" ini . basic assumption of Shellman and Morfitt [Ref. 5] that 

both the ginary parts of the modal function are linear along every 

edge of a itk is completely abandoned. For the evaluation of the A - t and 

B t coefficients, the "test for evanescence" conditions have been removed. A condition 
to determine whether to evaluate the coefficients from the ground level up or from 



7 



the top level down has been fomulated and incorporated into the program. This 
accomplishment leads to the relaxation of mode locating accuracy requirements 
which, combined with the improved precision of the revised program, makes the first 
order Newton-Raphson iteration unnecessary. The specific changes in the program 
and the resulting gains in speed, accuracy and execution stability are discussed in the 
following chapters. Suggestions to completely revise the mode search protocol to do 
without the "contour rectangles" and to look for the modes according to their range 
attenuation rates are also provided. 
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II. PROGRAM REVISIONS 



M-Layer is structured into three parts: setup, mode search and propagation 
factor evaluation. The main input is the modified refractive index values at specified 
heights so that a piecewise linear profile can be constructed. If the mode locations 
for the particular profile are available from a previous run of the program, they can 
also be included in the input and the mode search procedures will be bypassed. The 
various ranges and transmitter and receiver heights for which propagation factors are 
desired are also specified. The subroutine WVGSTDIN is called to input the 
information from an ASCII data file. The program then computes the constants to 
be used for mode search and propagation factor evaluation. The mode search is 
performed with the subroutine FNDMOD. The MODSUM subroutine is then 
invoked to first compute the and B { coefficients as explained in the Introduction, 
then compute the propagation factor and the propagation loss. The complete 
program structure is given in Figures 1 and 2. There are several other subroutines 
which are not included in these and other figures, such as DHORIZ for computing 
the horizon distance between a transmitter and a receiver for reference purpose; 
CHKMOD, a maintenance routine for removing zeroes from reported mode 
locations by older versions of the program; or A02H20, a routine to compute the 
atmospheric absorption coefficient due to oxygen and water vapor. They will not be 
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Figure 1 Original M-layer subroutines structure. 






Figure 2 Original M-layer subroutines structure (continued). 
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discussed as they do not contribute directly to the main purpose of this program of 
locating the modes and computing the propagation factor. 

The program structure has been altered as shown in Figures 3 and 4. Since the 
A t and coefficients have to be evaluated only once, they are now obtained through 
a call to the subroutine ABCOEF directly from the main program right after the 
modes are located. Several subroutines are dropped in this revision for various 
reasons: The subroutines NORME and NORMRE are eliminated because they are 
no longer needed due to the change in complex number representation; the 
subroutines NOMSHX, FDFDTX and DXDETR are not used because the modes 
are now located with adequate precision without further iteration; the subroutine 
ADDX is not listed separately because it is called only once and has been reduced 
to only a few lines which are placed where the subroutine is called in the original 
program. On the other hand, changes in the mode search algorithm require the 
addition of two new subroutines: SURFO is a modified and simpler version of SURF; 
ROOTS replaces QUAD. Due to the change in complex number representation, all 
subroutines listed below FNDMOD and MODSUM have been revised, including 
their input/output lists. But except for SURFO and ROOTS, the utilities of these 
subroutines are the same as those of the original ones. Descriptions of these 
subroutines can be found in the report by Yeoh [Ref. 4]. 

The most significant changes have been made in XCADD, XCDAIT and 
XCDAIG for adopting the complex exponent representation and improving 
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Figure 3 New M-layer subroutines structure. 
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Figure 4 New M-layer subroutines structure (continued). 
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computation speed and accuracy; in FZEROX, FINDFX, ROOTS and SURFO for 
stabilizing and simplifying the mode search algorithm; and in ABCOEF for 
implementing the criteria to determine the reliable manner for evaluating the^f,- and 
Bj coefficients. These changes are discussed in the sections below. The source code 
listings of the completely new subroutines XCADD and ROOTS and the significantly 
revised subroutines FZEROX and ABCOEF, which are compiled with Microsoft 
FORTRAN version 5.00, are attached as Appendices A through D. Validation of 
the revised program has been carried out at 9.6 GHz for all the 21 profiles listed in 
Yeoh [Ref. 4]. 

A. ADDITION SUBROUTINE 

XCADD is the subroutine implementing the addition of complex numbers 
under the representation by their exponents. Given the double precision complex 
numbers Zj and z 2 as the exponents of the addends, this subroutine returns the 
exponent of the sum. Since a double precision number has an accuracy' of 53 bits, 
if the real parts of z 1 and z 2 differ by more than 53 bits, the exponent of their sum 
will simply be the one of the greater real part. When cancellation becomes serious, 
the square root of the addends is factored out first. Then the four-term Taylor series 
expansions of the resulting reciprocals are summed. Since the leading term of the 
sum of the Taylor series is a good estimate of the sum of the reciprocals and the 
relative error of the four-term Taylor series sum is proportional to the fourth order 
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of this leading term, the threshold for invoking this interpolation procedure is set at 
the highest possible value of 2~ 14 allowed under double precision. Experimenting 
with this procedure shows that this interpolation improves accuracy as long as the 
threshold is set at a number between 2~ 24 and 2 ~ 14 . 

B. AIRY FUNCTION EVALUATION 

Similar to the original program, the evaluation of the Airy function adopted the 
algorithm prescribed by Schulten, et. al. [Ref. 6]. In the new program, changes are 
made to follow the advice of Schulten, et. al. concerning the region within which a 
Taylor series expansion, instead of the faster Gaussian quadrature, has to be used to 
achieve double precision accuracy. Other changes in implementing the algorithm are 
described below. 

1. XCDAIT 

Due to the similarity in their Taylor series coefficients, the Airy function 
and its derivative are evaluated within a single loop. The relative accuracy of the 
derivative of the Airy function is set at the double precision limit of 2 ~ 54 . 

2 . XCDAIG 

Six term Gaussian quadrature is used for evaluating the Airy function and 
its derivative outside the circle of radius 4.97 centered at (0.90, 2.80) on the complex 
plane. The use of four-term quadrature outside a radius of 15 from the origin 
suggested by Schulten, et. al. is not adopted. The six-term quadrature in this range 
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retains a higher accuracy while overall speed improvement by using both the four- 
term and the six-term quadrature appears to be minimal. 

C. MODE LOCATING 

As explained in the Introduction, the modes are located at the zeroes of the 
modal function. These zeroes are located on the upper complex q n plane. Here q n 
is the value of qj on the earth’s surface, which, according to Eq.(2) of Chapter I, is 
a linear function of p 2 . For a horizontally propagating mode, p/k is close to unity. 
The maximum range attenuation rate specified for the desired modes, which 
corresponds to a limit on the imaginary part of p, determines approximately the 
upper bound for the imaginary part of the q n complex plane to be searched for 
modes. The Shellman and Morffit mode search procedure first divides the search 
region horizontally into "contour rectangles" each of which spans 160 meshes along 
the real q n direction. A mesh is a square whose size is an adjustable parameter of 
the order 10~ 4 at 9.6 GHz for most of the cases considered herein. This parameter 
is determined by the frequency and the slope of the modified index of reflection in 
the lowest layer of the profile. The search commences at the top left corner of the 
"contour rectangle" whose left edge has a real coordinate value close to the 
difference of the real parts of the q n values, with the minimum modified index of 
refraction and the index near the surface substituted into Eq.(2) of Chapter I. After 
the search over the initial rectangle is completed, the program moves to search the 
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next rectangle until a specified maximum number of modes are found or a specified 
number of "contour rectangles" have been searched. 

The search for zeroes makes use of the fact that a real function changes sign 
when it crosses a simple zero. Since a zero of a complex valued function F(q) is 
where both its real part and imaginary part vanish, a necessary condition for a point 
to be a zero is that it is on the intersection of two curves defined by Im{F(q)} = 
0 and Re{F(q)} = 0. The program searches around a "contour rectangle" for a sign 
change in Im{F(q)} across an edge of a mesh bordering the side of the "contour 
rectangle" to determine that a line of Im{F(q)} = 0 has been encountered. The 
search then follows this line into the meshes within the "contour rectangle", checking 
each mesh to see if a curve Re{F(q)} = 0 enters the mesh under investigation. All 
these steps make use only of the assumption that the zeroes of the modal function 
are simple. Once both the curve Im{F(q)} = 0 and the curve Re{F(q)} = 0 are 
determined to be present within a mesh, the location of their possible interception 
is estimated. An algorithm for this estimate is required. 

Shellman and Morffit [Ref. 5] introduced a further assumption that the 
functions Re{F(q)} and Im{F(q)} are both linear along the edges of a mesh. Based 
on this assumption, they try to estimate the locations where the curve Im{F(q)} = 
0 enters and leaves a mesh square, and the location of q,^ if a curve Re{F(q)} = 0 
also enters the same mesh. It is obvious that information about the locations where 
the curves enter and leave the mesh square is not essential. Furthermore, in the 18 
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m duct height case, the scheme causes the search path to loop around four 
contiguous meshes until the search is broken up by the limit on the number of 
meshes to be investigated. Replacing their technique requires major changes in the 
subroutines involved. A new subroutine ROOTS is provided to estimate the location 
of the intersection of the curves Im{F(q)} = 0 and Re{F(q)} = 0. These changes 
eliminate the looping problem. 

Another problem is encountered in the 40 m duct height case when a large 
number of zeroes are found in the lower half complex q n plane. These zeroes 
appear to belong to the reflection coefficient on the wrong sheet of the branch cut 
and are not waveguide modes. This happens because the search region has been 
extended below the real q 11 axis to avoid the singularity in SURF. The problem with 
this singularity should have been solved within SURF, especially because it occurs 
only when the derivative of the subroutine output variable gamma with respect to q^ 
is computed. Since this derivative is not needed during mode search, the extension 
of the search region to the negative q n plane is unnecessary. A simplified routine, 
SURFO, is introduced which is exactly the same as SURF except that it does not 
evaluate the derivative of gamma. By using this subroutine instead of SURF, the 
search path in the revised program does not avoid the real and the imaginary axes. 

1. FNDMOD 

The search region is limited to the upper half q n plane. All the modes 
found are ordered according to their range attenuation rates before those numbered 
beyond the maximum modes allowed are abandoned. 
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2. FZEROX 



Since the curve Im{F(q)} = 0 enters into a mesh square through an edge, 
the values of Im{F(q)} must change sign over the end points of either one or all 
three other edges. When there is only one other edge across which Im{F(q)} 
changes sign at its end points, it is the edge across which the curve Im{F(q)} = 0 
exits the mesh square. Ambiguity arises when all edges indicate a change of sign at 
their end points. When this occurs, a "right turn rule" is adopted which assumes that 
the curve exits the edge to the right of the one along which it enters the mesh 
square. Such a rule avoids the retracing of the search path when the mesh square 
is revisited as entering this same mesh square from the left side of an edge after 
exiting from its right side requires a crossing of the Im{F(q)} = 0 curve, which is 
prohibited under the simple zero assumption. On the other hand, the actual curve 
may have turned left and then returns to this mesh square, i.e., following a "left turn 
rule." Under such a scenario, this wrong choice would have left a segment of the 
curve not searched. This difficulty has not been observed during testing. In fact, the 
ambiguous situation seldom occurs. Note also that, as remarked above, two lines of 
Im{F(q)} = 0 do not cross each other unless a higher order zero is present. Hence, 
only a "right turn rule" or a "left turn rule" for the curve to exit the mesh is allowed. 
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Exiting the opposite edge demands a pair of crossing Im{F(q)} = 0 curves within the 
mesh square. This violates the assumption that all zeroes are simple. Also note that, 
the possibility of vanishing Re{F(q)} or Im{F(q)} values at the corners of a mesh 
square is eliminated through a small adjustment in FINDFX. 

3. FINDFX 

Both the vertical shift away from the real q n axis and the horizontal 
offset away from the imaginary axis are unnecessary and have been removed from 
this routine. Furthermore, as a result of converting to the complex exponent 
representation, the sine and cosine of the argument of the modal functions are 
examined for sign changes in FZEROX. This is implemented in FINDFX by 
including the cosine and sine values of the argument of the modal function in the 
output list. To avoid the indeterminate case when either the real or the imaginary 
part of the modal function becomes zero at any corner of a mesh square, the 
argument for computing the cosine and sine values is increased by 2 -53 when this 
occurs. This is equivalent to a consistent small distortion of the particular corner of 
the mesh square. This will not cause any error in locating the zero because FINDFX 
still returns separately the unmodified exponent of the value of the modal function. 

4. ROOTS 

Assuming that the modal function is analytic within the mesh, this 
subroutine utilizes the values of the modal function at the four corners of the mesh 
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square to determine the Taylor series expansion coefficients of the modal function 
to the third order. The roots of this cubic polynomial are then located using 
Cardan’s solution by radicals. If the higher order coefficients fall below machine 
resolution for a root within the mesh square, these coefficients are regarded as zero 
and the order of the polynomial is reduced and can be solved more expediently. If 
the function is determined to be constant over the mesh square, the center of the 
square is taken as the root location. 

D. EVALUATING A- AND B ( 

As discussed in the Introduction, the A i and B i coefficients can be evaluated 
either from the top level down or from the lowest level up. These two procedures 
are simply called "integration down" and "integration up", respectively, in the original 
documentation [Ref. 4], The location of a mode has been called an eigenvalue. 
That the results of "integration down" and "integration up" agree is a manifestation 
that the eigenvalue is located accurately. 

The subroutine ABCOEF evaluates the coefficients A i and B for each mode. 
If the range attenuation rate for a mode is greater than 0.1 dB/km, the coefficients 
are evaluated from the lowest layer up. Otherwise, it is evaluated from the top layer 
down. It is obvious that such a rule must be implemented because the results of 
"integration up" and "integration down" do not agree for many modes. Efforts are 
made to determine the cause of this discrepancy and to devise a means to resolve it. 



20 



Investigation reveals that inadequate precision in the location of the modes is 
one source of the problem. Since the B i coefficients depend on the A i coefficients, 
while the A i coefficients are obtained directly, only the A i coefficients need to be 
examined. Th tA { coefficients of the six modes of lowest range attenuation rates for 
all 21 profiles except the one without evaporation duct are computed using 
eigenvalues of different accuracy controlled by the first order Newton-Raphson 
iteration method. Table 1 shows th eA- coefficient computed with the new program. 
They are arranged from the top layer down. In the i-th layer, the A i coefficient 
computed by "integration downward" depends only on A i+] in the layer above while 
that computed by "integration upward" depends only on A^j in the layer below. 
Hence in each layer, the coefficient obtained by "integration downward" is listed 
above that obtained by "integration upward". There are five sets of A i values listed, 
with the magnitudes given in powers of 10, and the phase given as a multiple of n. 
They are obtained from eigenvalues of decreasing accuracy -the one used to compute 
the left most column being the most accurate. The first set is computed using an 
eigenvalue having a relative accuracy of 2 -40 . The second set uses an eigenvalue with 
a relative accuracy of 2 -36 . The relative accuracy of the eigenvalue for the third set 
is 2 . For the fourth set, the first order Newton-Raphson iteration of the mode 

location is set at an absolute accuracy of 0.03 of the mesh size, same as that specified 
in the original program. The eigenvalue for the right most set is the mode location 
estimated by ROOTS without modification by the Newton-Raphson iteration. It is 
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TABLE 1. IMPROVING Aj ACCURACY WITH EIGENVALUE 18 M DUCT 



mode 


q- eigenvalue: 


. 18885743251768030+00 


.10806787448105980- 


■01 






eigenvalue difference: 


.000+00 


.00D+00 


-.490-13 - 


.660-15 


. 12D-1 1 - 


.120-10 


. 15D-06 


.600-07 


layer # 


A i /down 


A i /down 


A i /down 


A i /down 


A i /down 


layer # 


Ai/up 




Ai/up 


Ai/up 


Ai/up 


Ai/up 


18 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


18 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


17 


*.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


17 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


16 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


16 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


15 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


15 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


14 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


14 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


13 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


13 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


12 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3843 


.5659 


12 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


11 


2.2002 - 


.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.1909 


-.8068 


11 


2.2002 - 


.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


10 


5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4654 


.2423 


-4.1810 


-.2161 


10 


5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4647 


.2423 


9 


3.6974 - 


.6979 


-3.6974 


-.6979 


-3.6974 


-.6980 


-3.6783 


-.7012 


-6.4611 


-.2121 


9 


3.6978 - 


.6978 


-3.6978 


-.6978 


-3.6978 


-.6978 


-3.6978 


-.6978 


-3.6977 


-.6978 


8 


.3459 - 


.7982 


.3459 


-.7982 


.3460 


-.7982 


.3482 


-.7926 


-1.9078 


-.9148 


8 


.3459 - 


.7983 


.3459 


-.7983 


.3459 


-.7983 


.3459 


-.7983 


.3459 


-.7983 


7 


.4098 


.8794 


.4098 


.8794 


.4098 


.8794 


.4136 


.8836 


-1.0899 


.5364 


7 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


6 


.3480 


.8161 


.3480 


.8161 


.3480 


.8161 


.3526 


.8205 


-.5879 


.4005 


6 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


5 


.2923 


.8304 


.2923 


.8304 


.2923 


.8304 


.2972 


.8358 


-.3490 


.3749 


5 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


4 


.2359 


.8619 


.2359 


.8619 


.2360 


.8619 


.2408 


.8690 


-.2058 


.3731 


4 


.2358 


.8618 


.2358 


.8618 


.2358 


.8618 


.2358 


.8618 


.2358 


.8618 


3 


.1831 


.8910 


.1831 


.8910 


.1832 


.8910 


.1878 


.9003 


-.1250 


.3753 


3 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


2 


.1300 


.9149 


.1300 


.9149 


.1301 


.9149 


.1342 


.9275 


-.0734 


.3750 


2 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


1 


.0586 


.9335 


.0586 


.9335 


.0588 


.9335 


.0618 


.9545 


-.0318 


.3670 


1 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 
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clear that, for this mode, the difference between these two methods of computing the 
coefficients becomes negligible as the accuracy in mode location increases. For 
example, in the 8-th layer, the magnitude of A i computed by integrating downward 
changes from -1.9078 to 0.3482, to 0.3460 to 0.3459, which agrees with the result 
computed by integrating upward. The phase follows the same trend to an agreement 
within 0.001 7T. Table 2 shows a similar set of output, but the coefficients fail to agree 
even when the relative accuracy is increased to 2 _4 °. Note that the actual difference 
in both the real part and the imaginary part of the two most accurate eigenvalues is 
about 2 -48 . Double precision accuracy appears to be insufficient for the coefficients 
computed with these two methods to agree for all modes. Some interesting features 
can be observed in both tables, which are present in all 120 sets of values computed. 
When disagreement is present in one set of A i coefficients, such as those in either 
Table 1 or Table 2, the change toward smaller differences with improving eigenvalue 
accuracy occurs mainly in one way of computation, but not both. For example, in 
Table 1, the values of "integration downward" improve with better eigenvalue 
accuracy, while those computed by "integrating upward" change little. In Table 2, the 
results of "integration downward" are the ones that are holding steady as the accuracy 
in eigenvalue improves. Furthermore, when disagreement occurs, the layer in which 
the A i coefficient has the smallest magnitude, i.e., the one having the most negative 
power of 10, divides the table into two parts. The results of two different ways of 
computation agree in the layers above this one if they disagree in those below it, and 
vise versa. No explanation will be attempted. Instead, practical rules are drawn up 
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TABLE 2. IMPROVING Aj ACCURACY WITH EIGENVALUE 36 M DUCT 



mode 


3 q- eigenvalue 


.31 480001 64781 392D+00 


.14796229405720070- 


•02 






eigenvalue difference: 


.380-14 


.360-14 


.380-14 , 


.360-14 


-.160-09 - 


.220-09 


- .53D-07 . 


.280-07 


layer # A i /down 


A i /down 


A i /down 


A i /down 


A i /down 


layer # Ai/up 




Ai/up 


Ai/up 


Ai/up 


Ai/up 


27 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


27 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


26 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


26 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


25 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


25 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8852 


.3913 


24 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


24 


-7.4914 


.6081 


-7.4914 


.6081 


-7.4914 


.6081 


-7.4914 


.6081 


-7.4914 


.6081 


23 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


23 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


22 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


22 


-23.5827 - 


.9407 


-23.5827 


-.9406 


-23.5827 


-.9406 


-23.5827 


-.9406 


-23.5827 


-.9406 


21 


-35.2395 - 


.2502 


-35.2395 


-.2502 


-35.2395 


-.2502 


-35.2395 


-.2502 


-35.2396 


-.2501 


21 


-44.4517 - 


.8199 


-45.8691 


.8599 


-45.8691 


.8599 


-47.4590 


-.1252 


-47.4594 


-.1251 


20 - 


131.3304 - 


.9570 


-131.3304 


-.9570 


-131.3304 


-.9570 


-131.3304 


-.9570 


-131.3307 


-.9569 


20 - 


129.0146 - 


.2961 


-127.6070 


.0248 


-127.6070 


.0248 


-122.9124 


-.9081 


-120.9305 


.8279 


19 


-25.6088 - 


.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


19 


-25.6090 - 


.9228 


-25.6184 


-.9241 


-25.6184 


-.9241 


-22.5644 


-.8054 


-20.2166 


.7391 


18 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


18 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


-13.0618 


.7675 


-10.8148 


.3440 


17 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


17 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0308 


.4179 


-6.3129 


.1800 


16 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


16 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3116 


.2970 


15 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


15 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3127 


.2629 


14 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


14 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5668 


.2415 


13 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


13 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


12 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


12 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


11 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


11 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


10 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


10 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 
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to take advantage of these facts. In Table 1, the process of "integration upward" goes 
through the troublesome 10-th layer and produces results which agree with the results 
of downward integration before the downward process goes through the 10-th layer. 
On the other hand, the downward integration is tripped up going across the 10-th 
layer and produces results which fail to agree with the results from upward 
integration. It is clear that the results from upward integration are the correct ones. 
This conclusion is further supported by the fact that improving the accuracy of the 
eigenvalue does not change significantly the results of upward integration. Similar 
argument leads to the conclusion that in Table 2, the results of downward integration 
are the correct values. 

It can be concluded from the above observations that one of the methods of 
computing the A i coefficients converges to the correct value much faster than the 
other. It is also found that this method of faster convergence is always able to arrive 
at the correct values for A i for all the cases under investigation. 

Table 3 lists the statistics of the method of integration which yields the correct 
A i coefficients for each of the 120 modes investigated. The differences in magnitudes 
and phases in the lowest layer and in the layer below the highest are also listed. 
Since for most of the cases, when disagreement in A i values occurs, the correct 
integration is upward -this is used as the default. To decide that downward 
integration should be utilized, the following steps are taken: The first A i value of 
downward integration is computed and compared to the value from upward 
integration. If the magnitudes in dB disagree by less than 0.02 dB, their phases wall 
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be checked. If the phases differ by less than 10 -3 7r, the agreement is deemed 
acceptable and the A i and B i coefficients computed from the lowest layer up are 
used. Otherwise, the coefficients are re-evaluated again from the highest layer down. 

Once the correct method of evaluating the A i and B- t coefficients is used, the 
accuracy of the mode location becomes less critical. For all the cases investigated, 
the A i coefficients obtained from mode locations estimated with or without the 
Newton-Raphson first order iteration differ only by 0.06 dB in magnitude and 
0.0013 7r in phase at most. In fact, few cases show differences more than 0.002 dB 
and 0.0001 ir. The Newton-Raphson iteration is not needed. Hence the subroutines 
NOMSHX, FDFDTX and DXDETR are removed. 
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TABLE 3. STATISTICS FOR EVALUATING \ COEFFICIENT 



Duct 

height 


Mode # 


Evaluating Method 


> 

> 

i 


A arg (A,)/ n 

f 1 1 


Layer 


Layer 


up 


down 


bottom 


top-1 


bottom 


top-1 


02 


1 


X 












2 


X 












3 


X 












3 


X 




0.172 




0.093 




5 


X 












6 


X 




8362 




13234 




04 


1 


X 












2 


X 












3 


X 




0.008 




0.0002 




3 


X 




1.030 




1.8717 






5 


X 




7.814 




13948 




6 


X 




0.002 




0.0001 




06 


1 


X 












2 


X 




0.002 




0.0004 




3 


X 




0.522 




0.0158 




4 


X 












5 


X 




13278 




0.4377 




6 


X 




0.002 




0.0001 




08 


1 


X 












2 


X 




0.002 








3 


X 




0.002 




0.0001 


0.0001 


4 


X 




0.016 




0.0026 




5 


X 




4.066 




0.6355 




6 


X 




3.978 




0.6186 
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TABLE 3. CONTINUED 1 



Duct 


Mode# j 


Evaluating Method 


A 141 


m 


bargiAJ/n 


height 










! 


















Layer 


Layer 






up 


down 


bottom 


top-1 


bottom 


top-1 




1 


X 














2 


X 








0.0002 






3 


X 








0.0001 




10 


4 


X 




0.04 




0.0008 






5 


X 




0.206 




0.0402 


0.0001 




6 


X 




0.002 




0.0*08 






1 


X 














2 


X 




0.006 




0 . 000 * 






3 


’ X 




0.004 








12 


4 


X 




1.808 




0.5661 






5 


X 




1.732 




0.5429 






TJ 


X 




1.472 




0.0414 






1 


X 














2 


X 




0.102 




O.IOOS 






3 


X 




0.178 




0.0052 




14 


4 


X 




0.024 




0.0005 






5 


X 




0.004 




0.0001 






TT 


X 




0.85 




0.4711 






1 


4 














2 


4 




0.006 




0.0002 






3 


X 




0.004 








16 


4 


X 




0.006 




0.0001 






5 


X 




0.002 




0.0001 






6 


X 




0.004 




0.0077 


* 
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TABLE 3. CONTINUED 2. 











AKI 


m 


bargiAJIn 


Duct 


Mode # 


Evaluating Method 










height 




































Layer 


Layer 






up 


down 


bottom 


top-1 


bottom 


top-1 




1 




X 




0.008 




0.0001 




2 


4 




0.002 




0.0001 






3 


X 








0.0001 




18 


4 


X 














» 


X 




0.016 




0.0003 






6 


4 




0.002 










1 




4 




0.078 




0.0164 




2 


X 














3 


X 




0.002 




0.0001 




20 


4 


X 








0.0008 






5 


4 




0.16 




0.0195 






6 


4 




0.002 




*.0001 






1 




X 




8.708 




0339 




2 


X 














3 


X 




0.004 








22 




4 




0.016 










6 


X 




0.002 




0.0001 






6 


X 




031 




0.0117 






1 


X 














2 




X 




0.868 




03842 




3 


X 




0.006 




0.0009 




24 


4 


X 




0.002 




0.0001 






5 


X 




0.026 




0.0009 






6 


X 




0.008 




0.0001 





29 



TABLE 3. CONTINUED 3 



Duct 


Mode # 


Evaluating Method 


AM,I 


1 


A arg (X,.)/ it 


height 














1 






















Layer 


Layer 






up 


down 


bottom 


top-1 


bottom 


top-1 




1 


X 




0.002 


0.002 


0.0002 


0.0001 




2 




X 




4308 




0.121 




3 


X 




0.006 








26 


4 


X 




0.002 




0.0001 






5 


X 








0.0001 






6 


6 




0.034 




0.0039 






1 




X 




0.028 




0.0014 




2 




X 




4.806 




0.0728 




3 


rX 












28 


3 


X 














5 


X 




0.004 


0.002 


0.0002 






6 


X 




0.004 




0.0019 






1 




X 




1.562 




0.0165 




2 


6 














3 




X 




0.718 




03455 


30 


X 


6 














5 


X 




0.004 










6 


X 




0.724 




0.0522 






1 




X 




3.194 




0.1648 




2 


X 




0.002 










3 




X 




13.12 




0.1026 


32 


4 


X 




0.002 










5 


X 




0382 




0.0099 






_6 ... 


X 




0.002 




0.0001 
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TABLE 3. CONTINUED 4. 











AKI 


m 


A arg(A t )l it 


Duct 


Mode # 


Evaluating Method 










height 




































Layer 


Layer 






up 


down 


bottom 


top-1 


bottom 


top-1 




1 


X 




0.002 


0.002 








2 




X 




13.456 




0.0311 




3 




X 




1.014 




02347 


34 


4 


X 














5 


X 




0.03 




0.0006 






6 


X 




0.014 




0.0006 






1 




X 






0.0001 


0.0014 




2 




X 




1.686 




02224 




4 




X 




4.724 




0.0919 


36 


4 


X 














6 


X 




0.006 




0.0001 






6 


X 




0.02 




D.0I01 






1 




X 




0.996 




0.0115 




2 




X 




4.974 




0.0152 




3 


X 












38 


4 




X 




5.052 




0.0417 




5 


X 








9.00I1 






6 


X 




0.002 










1 


X 




0.002 


0.002 








2 




X 




3.85 




0.1226 




3 




X 




3.568 




0.1555 


40 


4 




X 




3.448 




0.1678 




5 


X 








0.0001 






- 6 


X 













31 



III. CONCLUSIONS AND RECOMMENDATIONS 



A. Performance 

This revision of M-Layer converts the extended complex number representation 
of an exponentially large or small number into the direct representation by its 
complex exponent. The accuracy of the computation has been improved in two ways: 
First, an interpolation algorithm has been devised when severe cancellation of the 
addends is detected. Secondly, accuracy for the evaluation of the Airy function has 
been improved, not just by summing the Taylor series to double precision resolution 
and by adopting six-term Gaussian quadrature, but also by expanding the region 
within which the more expedient Gaussian quadrature is excluded in favor of the 
more accurate, but time-consuming, Taylor series summation. The improvement in 
accuracy is most easily seen from Table 1. 

As discussed in the Introduction, evaluating the Aj and B - t coefficients either 
from the lowest layer up (integration up), or from the top layer down (integration 
down), must result in the same values. This property provides a consistency check 
for the accuracy of the computation. For the six modes of lowest range attenuation 
rates of the 20 profiles of different duct heights, Table 1 lists the maximum 
difference for each mode which shows a discrepancy between these two methods of 
evaluating the A t - coefficients. For each profile, the maximum value in magnitude 
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TABLE 1. MAXIMUM DIFFERENCE IN \ COEFFICIENT BETWEEN 
INTEGRATION UP AND DOWN 



Duct 

height 


Mode 

# 


Difference in A coefficient 


Magnitude difference in 
(dB) 


Phase difference over 

0.1 7 r 


old 


new 


old 


new 


02 


6 


5.22 




Yes 




6 


61.16 




Yes 




04 


6 


22.46 


2.3 






5 


106.9 




Yes 




06 


4 


8.62 




Yes 




5 


32.36 








02 


5 


77.84 




Yes 




6 


44.9 




Yes 




14 


6 






Yes 




12 


4 


69.38 




Yes 




5 


46.32 




Yes 




6 


7.46 




Yes 




14 


6 


30.6 




Yes 




22 


6 


8.64 




Yes 




24 


2 


80.48 




Yes 




26 


2 


110.68 




Yes 




28 


2 


1§6.9 


67.68 


Yes 


Yes 


30 


3 


173.28 


143.42 


Yes 


Yes 


32 


1 


11.38 




Yes 




3 


525.04 


188.04 


Yes 


Yes 
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TABLE 1. CONTINUED 



Duct 

height 


Mode 

# 


Difference in A coefficient 


Magnitude difference in 
(dB) 


Phase difference over 

0.1 7 r 


old 


new 


old 


new 


34 


2 


37.98 




Yes 




4 


715.7 


209.94 


Yes 


Yes 


34 


2 


112.74 




Yes 




4 


957.92 


231.68 


Yes 


Yes 


38 


2 


107.44 


52.26 


Yes 


Yes 


4 


1249 


255.8 


Yes 


Yes 


34 


3 


167 


112.72 


Yes 


Yes 


4 


823.56 


258.18 


Yes 


Yes 




Magnitude difference within 2dB are not listed. 



difference in dB among all the layers is listed if it is greater than 2. If the phases of 
the coefficients deviate more than O.Itt in any layer, that particular mode is also 
singled out. The location of the mode of the revised program is within a relative 
accuracy of 2 -40 achieved through first order Newton-Raphson iteration. Even 
though discrepancies still exist when the duct is 28 meters or higher, it is clear that 
the revised program computes more accurately than the original one. 

For the cases where the two methods of evaluating the A i and B i coefficients 
disagree, it has been observed that one of the methods always leads to A i values 
which are little changed when the accuracy in mode location is varied, while the 
other method produces A i values which shift toward the results of the other method 
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as the accuracy of mode location improves. Based on this observation, a consistency 
check is implemented into the program to identify the method which converges 
better. For the 120 cases investigated, when this method of faster convergence is 
used, the ^coefficients obtained from mode locations estimated with or without the 
Newton-Raphson first order iteration differ only by 0.06 dB in magnitude and 
0.0013 n in phase at most. In fact, few cases show differences more than 0.002 dB 
and 0.0001 n. This allowed the Newton-Raphson iteration to be removed in this 
revision. 

Table 2 compares the performance between the original and the revised 
programs. The time spent to find the modes has been reduced by an average of 
22.58%. The revised program can always produce the modes found by the original 
program. Moreover, the mode search is stable for the new program: the time it 
requires to search for the modes is about the same for similar profiles. The sudden 
jumps in mode search time for the 24 m and the 40 m cases, which indicate troubles 
during the search, no longer happen. 

With the proper method of evaluating the A { and B i coefficients determined by 
the consistency check, the output of the revised program differs from the original 
program in some cases. The most serious deviation has been observed for the 38 m 
duct height case as shown in Tables 3 and 4. For example, at a range of 36.5 km 
with the transmitter at a height of 25 m and the receiver at 10 m, the coherent path 
loss is 175.93 dB from the original program, and is 167.90 dB from the revised 
program. 
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TABLE 2 OVERALL MODE SEARCH PERFORMANCE COMPARISON 



DUCT 

HEIGHT 

(meter) 


ORIGINAL 


REVISED 


Time 

Improvement 


Time 


Modes 


Time 


Modes 


40 


0:00:37 


3 


0:00:35 


3 


5.40% 


02 


0:32:14 


9 


0:31:55 


9 


0.98% 


04 


1:14:12 


25 


1:05:04 


25 


12.31% 


05 


2:10:18 


53 


1:56:50 


53 


10.33% 


08 


0:35:58 


39 


0:29:25 


39 


18.21% 


10 


0:53:24 


59 


0:48:32 


61 


9.11% 


12 


1:09:40 


86 


1:01:44 


89 


11.39% 


14 


1:20:42 


94 


1:11:13 


97 


11.75% 


15 


1:54:35 


95 


1:46:47 


97 


31.82% 


40 


1:45:09 


109 


1:27:15 


104 


17.02% 


24 


1:46:19 


10* 


1:34:20 


105 


11.27% 


22 


1:52:54 


105 


1:35:18 


106 


15.59% 


24 


3:42:59 


109 


1:46:47 


107 


52.11% 


25 


2:07:42 


109 


1:43:55 


106 


18.62% 


28 


2:00:05 


107 


1:44:59 


109 


12.57% 


30 


1:59:59 


107 


1:46:47 


106 


11.39% 


32 


1:55:29 


108 


1:42:58 


110 


10.84% 


34 


2:29:57 


109 


2:15:58 


111 


9.32% 


36 


2:31:40 


109 


2:17:20 


112 


9.45% 


30 


2:38:44 


110 


2:18:09 


111 


12.97% 


40 


5:41:17 


95 


2:39:39 


111 


53.22% 


Total 


40:23:54 




31:16:22 




22.58% 



36 



TABLE 3. ORIGINAL PROGRAM 38 M DUCT OUTPUT 



frequency = 


9600. 000C 


mhz 










range 


zt 


zr 


coherent 


incoherent 


coherent 


incoherent 


horizon 


(km) 


(m) 


(m) 


mode sum 
(dB) 


mode sum 
(dB) 


path Loss 
(dB) 


path Loss 
(dB) 


(km) 


27.3 


25.0 


4.0 


-15.30 


-15.62 


156.10 


156.43 


28.9 


27.3 


25.0 


6.0 


.62 


-2.35 


140.18 


143.16 


30.7 


27.3 


25.0 


8.0 


-1.11 


-4.21 


141.92 


145.01 


32.3 


27.3 


25.0 


10.0 


-27.26 


-12.66 


168.06 


153.46 


33.6 


36.5 


25.0 


4.0 


-16.94 


-16.62 


160.28 


159.96 


28.9 


36.5 


25.0 


6.0 


-.73 


-2.05 


144.07 


145.39 


30.7 


36.5 


25.0 


8.0 


-2.21 


-3.72 


145.55 


147.06 


32.3 


36.5 


25.0 


10.0 


-32.59 


-14.29 


175.93 


157.64 


33.6 


45.8 


25.0 


4.0 


-19.89 


-16.96 


165.20 


162.26 


28.9 


45.8 


25.0 


6.0 


-2.81 


-1.89 


148.11 


147.19 


30.7 


45.8 


25.0 


8.0 


-4.11 


-3.43 


149.41 


148.74 


32.3 


45.8 


25.0 


10.0 


-28.57 


-15.22 


173.88 


160.52 


33.6 



TABLE 4. REVISED PROGRAM 38 M DUCT OUTPUT 



frequency = 


9600. 000C 


) mhz 










range 


zt 


zr 


coherent 


incoherent 


coherent 


incoherent 


hori zon 


(km) 


(m) 


(m) 


mode sun 
(dB) 


mode sum 
(dB) 


path loss 
(dB) 


path loss 
(dB) 


(km) 


27.3 


25.0 


4.0 


-14.38 


-15.66 


155.18 


156.47 


28.9 


27.3 


25.0 


6.0 


.42 


-2.37 


140.39 


143.18 


30.7 


27.3 


25.0 


8.0 


-1.52 


-4.21 


142.33 


145.02 


32.3 


27.3 


25.0 


10.0 


-21.20 


-12.51 


162.01 


153.31 


33.6 


36.5 


25.0 


4.0 


-17.32 


-16.60 


160.66 


159.94 


28.9 


36.5 


25.0 


6.0 


-.48 


-2.08 


143.82 


145.42 


30.7 


36.5 


25.0 


8.0 


-1.62 


-3.73 


144.96 


147.07 


32.3 


36.5 


25.0 


10.0 


-24.56 


-14.04 


167.90 


157.38 


33.6 


45.8 


25.0 


4.0 


-20.26 


-16.93 


165.57 


162.23 


28.9 


45.8 


25.0 


6.0 


-3.14 


-1.93 


148.44 


147.23 


30.7 


45.8 


25.0 


8.0 


-4.62 


-3.46 


149.92 


148.76 


32.3 


45.8 


25.0 


10.0 


-25.40 


-14.90 


170.71 


160.21 


33.6 
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B. Recommendations 



The mode search protocol of this program needs to be revised. Since the 
search is limited by the number of modes to be found and the maximum range 
attenuation rate accepted, it is more logical to begin with locating the mode of the 
lowest attenuation, and then proceed to look for the next one in the order of 
increasing attenuation rate. Furthermore, there appears to be only a single ‘phase 
line’ of vanishing real part of the modal function on which all the modes are located. 
This line extends from lower to higher range attenuation rates. The partition of the 
search region into rectangles, as has been done in this program, tends to cut the 
"phase line" into segments before the program starts to search for the end points of 
these segments and then follow the segments in different directions. It is clear that 
a better way for mode search is to find the lower end of the single "phase line" then 
follow it to the other end. 
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APPENDIX A: SUBROUTINE XCADD 



This Appendix lists the addition subroutine XCADD which returns the complex 
exponent of the sum when the complex exponents of the addends are given. This is 
a complete re-write of the original subroutine of the same name. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 



subroutine xcadd(zx, zlx, z2x) 



c Given zlx and z2x, this subroutine adds the two complex numbers 
c z1=exp(z1x) and z2=exp(z2x) for z=exp(zx) and returns zx. 

c 

c inputs... 

c z1x=complex exponent of the complex number zl 

c z2x=complex exponent of the complex number z2 

c 

c outputs 

c zx=complex exponent of the complex number z 

c 

c subroutines called... 

c 

£**★★★★**** 

implicit real*8 (a-h,o-z) 

complex* 16 zx, zlx, z2x,zt1x,zt2x,clogzh,dsum,czero,cerrx,cone,chpi 
parameter ( pi =3 . 1 4 1 592653589793238462643d0 , t wopi =2 .d0*pi , 

+ hpi =0.5d0*pi ,zero=0.d0,c16=1 .dO/6.dO, 

+ bit 14=1. dO/16384. d0,bit24=bit 14/1024. d0,ctol=bit 14, 

+ dpi =2259 . d0/4294967296 . d0/4294967296 . dO , hdp i =dpi /2 . dO , 

+ e2m54=-3. 74299477502370481 9d1 ,e2p27=-0 . 5d0*e2m54 , 

+ chp i = ( 0 . dO , 1 . 57079632679489661 923 1 32d0 ) , cone= ( 1 . dO , 0 . dO ) , 

+ czero=(0.d0,0.d0) ,cerrx=( -3. 742994 77502370481 9d1 ,0.d0)) 

c cerrx=e2m54=-54*log(2)=exponent below machine accuracy 

dimension ztmp(2), stmp(2) 
equivalence (ztmp,clogzh),(stmp,dsum) 

£**★** 

c Replace the input variables with a local variable so that 
c equations in the form of y=x+y will not lead to confusion, 
c 

zt1x=z1x 

zt2x=z2x 

c 

clogzh=0.5d0*(zt1x-zt2x) 
dxh=ztmp(1 ) 

if(dxh .It. zero) then 
zx=zt2x 
dxh=-dxh 
else 

zx=zt1x 
end if 
£*★★★★★★★ 

c machine accuracy = 2**(-53) 

c 2**(27)=e**e2p2 7 
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46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 



if (dxh .ge. e2p27) then 
return 
else 

zx=0.5d0*(zt1x+zt2x) 
dsum=cdexp(clogzh) 
dsum=1 .dO/dsum+dsum 
if (cdabs(dsum) .gt. ctol) then 
zx=cdlog(dsum)+zx 
else 

Cancellation is serious. Imlclogzh] is close to pi/2 or -pi/2. 
yi=dnint(ztmp(2)/twopi )*2.d0 
ztmp(2)=ztmp(2)-pi*yi 
dyi =dpi*yi 

if (ztmp(2) .It. zero) then 
clogzh=-clogzh 
dyi=-dyi 
end if 

ztmp(2)=(ztmp(2)-hpi )-hdpi - dyi 

dsum=2 . dO* c l og zh * ( c one+ c 1 6*c l og zh * c l og zh ) 

if (dsum .eq. czero) then 

Note that a complete cancellation of two nonzero numbers of 
order one is considered to be as accurate as what is allowed 
by the machine and the algorithm. 
zx=cerrx+chpi+zx 
else 

dsum=cd log (dsum) 

if (stmp(l) .It. e2m54) stmp( 1 )=e2m54 
zx=dsum+chpi+zx 
end if 
end if 
return 
end if 
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APPENDIX B: SUBROUTINE FZEROX 



This Appendix includes the listing of the subroutine FZEROX which searches 
identifies the meshes which may contain modes within a contour rectangle. The 
Shellman-Morffit mode locating algorithm has been completely replaced. 
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1 subroutine fzerox(t left, tright, tbot, ttop, tmshO, zeros, ni ,nf ) 

2 £* * * * * 

3 c fzerox is a routine for finding the zeroes of a complex function, f , 

4 c which lie within a specified rectangular region of the 

5 c complex ql 1 plane, assuming that the function has only 

6 c simple zeroes over this rectangle. 

7 c 

8 c parameters specifying the search rectangle: 

9 c tleft - value of the real part of ql 1 at the left edge. 



10 


c 


tright- 


value of the real part of ql 1 at the right edge. 


11 


c 


tbot - 


value of the imaginary part of ql 1 at the bottom edge. 


12 


c 




(this is set to 0.) 


13 


c 


ttop - 


value of the imaginary part of ql 1 at the top edge. 


14 


c 


tmesh - 


set equal to about half the average spacing between 


15 


c 




zeroes within the rectangle. A smaller value may be used 


16 


c 




as a safety measure, but too small a value will result 


17 


c 




in excessively long run time. 


18 


c 


zeros - 


output list of (complex) values of ql 1 at which 


19 


c 




zeroes are found. 


20 


c 


nf-ni - 


the number of zeroes found 


21 


c 






22 


c 


subroutines 


. cal ledd- - 


23 


c 


f indfx 


24 


c 


roots 


25 


c 


nomshx 


26 


c* 


**** 




27 




impl ici t 


: double precision (a-h,o-z) 


28 




complex* 


16 f 1 0, f 01 , f 1 1 , fxnew, fxold, fxOO, fxlO, fxOI , f x 11, 


29 




+ 


czero,one, ci , sol , zeros 


30 




parameter(czero=(0.d0, O.dO) ,one=(1 .d0,0.d0),ci=(0.d0, 1 . dO ) ) 


31 


$include: 'mlaparm.inc' 



***** Begin listing of: mlaparm.inc 

1 c 

2 c include file to define the 

3 c maximum # of layers (mxlayr) 

4 c maximum # of modes (mxmode) 

5 c 

6 parameter (mx layr=35 ) 

7 parameter (mxmode=127) 

***** End listing of: mlaparm.inc 

32 dimension kedgel ( 100) , kedge2( 100) , kedge3(100) , kedge4( 100) , 

33 c + loc12r(mxmode) , loc12i (mxmode) , l oc 23 r (mxmode) , loc23i (mxmode) , 

34 c + l oc34 r ( mxmode ), l oc34i (mxmode), l oc41r (mxmode), l oc41 i (mxmode), 

35 + sol(3), theta(2) , zeros(2*mxmode+1 ) 

36 c 
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37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

73 

76 

77 

78 

79 

80 

81 



c 

common /tmccom/tmesh 

c***** 

c maxnsq - maximum number of mesh squares allowed on any one 
c phase line 

c maxnt - maximum number of times fzerox will reduce tmesh 
c 

maxnsq=3*max0( i nt ( ( t t op- t bot )/ tmshO ) , i nt ( ( t r i gh t - t l ef t )/ tmshO ) ) 
maxnt=2 

£*♦*** 

tmesh = tmshO 
ntime = 0 
go to 7 
c 

5 tmesh=tmesh/2.0d0 

ntime = ntime+1 
if (ntime .gt. maxnt) go to 97 
c 

7 continue 

★★ 

c calculate coordinates of rectangle edges in tmesh units 

c 

jit = idnint( tlef t/tmesh-0.5d0) 
jrt = idnint(tright/tmesh+0.5d0) 
jtop = idnint ( ttop/tmesh+1 .5d0) 
j bot = 0 
c 

c initialize parameters for starting search at upper left 

c corner of search rectangle 

c 

ki = jtop 
kr = jit 
kedge = 1 

call f indfx(kr # ki , fxnew, xnew,ynew) 

nre1=0 

nre2=0 

nre3=0 

nre4=0 

knot 12=0 

knot23=0 

knot34=0 

knot41=0 

nf=ni 

ni 1 =ni+1 
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82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 



go to 15 



£***★* 

10 continue 

i f (nrzl . It . 2) go to 15 
c wri te(16,2000) nrzl 

go to 5 
15 nrzl=0 

nrsqu = 0 

20 fxold=fxnew 
xold=xnew 
yold=ynew 

go to (21 ,26,31 ,36), kedge 

£*★*** 

c search along left edge of rectangle for changes in the 

c sign of imag(f) 

c 

21 continue 

i f ( k i .eq. jbot ) then 
kedge=2 
go to 26 
end i f 
ki = ki-1 

call f indfx(kr,ki , fxnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
if (nrel .eq.0) go to 23 
c 

c check if crossing point has been previously found 

c 

do 22 k=1,nre1 

i f(ki .eq.kedgel(k)) go to 20 

22 continue 
c 

c follow phase line through rectangular region 

c 

23 fx01=fxold 
fx01r=xold 
fxOI i=yotd 
f x00=f xnew 
fx00r=xnew 
fxOOi =ynew 
li = ki 

lr = jit 
go to 43 

c search along bottom edge of rectangle for changes in the 
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127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 
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158 

159 

160 

161 

162 
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164 

165 

166 

167 

168 

169 

170 

171 



c sign of imag(f) 
c 

26 continue 
if(kr.eq. jrt) then 

kedge=3 
go to 31 
end if 
kr = kr+1 

call f indfx(kr,ki,fxnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
i f (nre2.eq.O) go to 28 
c 

c check if crossing point has been previously found 

c 

do 27 k=1,nre2 
if (kr.eq.kedge2(k)) go to 20 

27 continue 
c 

c follow phase line through rectangular region 

c 

28 fx00=fxold 
fx00r=xold 
fx00i=yold 
fx10=fxnew 
fx10r=xnew 
fx10i=ynew 
li = jbot 
lr = kr-1 
go to 48 

c search along right edge of rectangle for sign changes in imag(f). 
c 

31 continue 

i f (ki .eq. jtop) then 
kedge=4 
go to 36 
end if 
ki = ki+1 

cal l f indfx(kr,ki ,fxnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
if (nre3.eq.O) go to 33 
c 

c check if crossing point has been previously found 

c 

do 32 k=1,nre3 



46 



172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 
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184 

185 

186 
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188 

189 
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191 
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193 
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195 
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197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 



if (ki .eq.kedge3(k) ) go to 20 

32 continue 
c 

c follow phase line through rectangular region 

c 

33 fx10=f xold 
fx10r=xold 
fx10i=yold 
fxl 1=fxnew 
fxl 1 r=xnew 
fx11i=ynew 
li = ki-1 
Ir = jrt-1 
go to 53 

£* ★★★★ 

c search along top edge of rectangle for sign changes in imag(f). 
c 

36 continue 

if (kr.eq. j l t) go to 80 
kr = kr-1 

call f indfx(kr,ki , fxnew, xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
i f (nre4.eq.O) go to 38 
c 

c check if crossing point has been previously found 

c 

do 37 k=1 , nre4 
if(kr.eq.kedge4(k)) go to 20 

37 continue 
c 

c follow phase line through rectangular region 

c 

38 fx11=fxold 
fxl 1r=xold 
f x 1 1 i =yold 
fx01=fxnew 
fxOI r=xnew 
f xOI i =ynew 
li = jtop-1 
lr = kr 

go to 58 

£***★* 

c enter mesh square from left side or exit rectangle at right edge. 
41 l r=l r+1 
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217 

218 

219 

220 

221 

222 

223 

224 

225 
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234 

235 

236 

237 

238 
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241 
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248 
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251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 



if (Ir .le. jrt-1) go to 42 
nre3=nre3+1 
kedge3(nre3)=l i + 1 
go to 10 

42 fx01 = fx1 1 
fx01r=fx1 1 r 
fxOI i=fx1 1 i 
fx00=fx10 
fx00r=fx10r 
fx00i=fx10i 

43 continue 

cell findfx(lr+1,li+1,fx11,fx11r,fx11i) 
call findfx(lr+1, li ,fx10,fx10r,fx10i ) 
c ******* 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei t=f xOI i*fx1 1 i 
edgeib=fx00i*fx10i 
if (edgeib .gt. O.dO) then 
c lm(f)=0 goes through the 01 to 10 line, 

if (edgeit .gt. O.dO) then 

c lm(f)=0 goes through the 10 to 11 edge (edge 1). 

lout=1 
else 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 
end if 
else 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

lout=4 

if (edgeit .It. O.dO) then 

c lm(f)=0 also runs through 01 to 11 and 10 to 11 edges, 

c Store crossing location and in/out information. 

knot34=knot34+1 
c loc34r(knot34)=lr 

c loc34i (knot34)=l i 

end if 
end if 
£* ★ ★ ★ ★ ★ ★ 

go to 60 

Q-kwkkk 

c enter mesh square from bottom side or exit rectangle at top edge. 
46 l i = l i+1 

if (li .le. jtop-1) go to 47 

nre4=nre4+1 

kedge4(nre4)=lr 
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263 
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267 
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301 
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go to 10 
f x00=f xOI 
fx00r=fx01r 
fxOOi =f xOI i 
f x10=f xl 1 
fx10r=fx1 1 r 
fx10i=fx1 1 i 
48 continue 

call findfx(lr,li+1,fx01,fx01r,fx01i) 
call f indfx(lr+1 , l i+1 , fxl 1 , f xl 1 r, fxl 1 i ) 
c ******* 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei l=fx00i*f xOI i 
edgei r=fx10i*fx1 1 i 
if (edgeir .gt. O.dO) then 
c lm(f)=0 goes through the 00 to 11 line, 

if (edgei l .gt. O.dO) then 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 
else 

c lm(f)=0 goes through the 00 to 01 edge (edge 3). 

lout=3 
end if 
else 

c lm(f)=0 goes through the 10 to 11 edge (edge 1) 

lout=1 

if (edgei l .It. O.dO) then 

c lm(f)=0 also runs through 00 to 01 and 01 to 11 edges, 

c Store crossing location and in/out information. 

knot41=knot41+1 
c loc41r(knot41)=lr 

c loc41 i (knot41 )=l i 

end if 
end if 

Q******* 

go to 60 

£***** 

c enter mesh square from right side or exit rectangle at left edge. 

51 lr=lr-1 

if (Ir .ge. jit) go to 52 
nre1=nre1+1 
kedgeKnrel ) = li 
go to 10 

52 fxl 1=fx01 
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307 

308 

309 
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335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 



fxl 1 r=fx01 r 
fxll i=fx01 i 
fx10=fx00 
fx10r=fx00r 
fx10i=fx00i 
53 continue 

call f indfx( Ir, l i + 1 , fxOI f fxOI r, fxOI i ) 

cal l f indfx( lr, l l ,fx00, fxOOr, fxOOi ) 

★ ★ ★ ★ ★ ★ 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei t=fx01 i*fx1 1 i 
edgeib=fx00i*fx10i 
if (edgeit .gt. O.dO) then 
c lm(f)=0 goes through the 01 to 10 line, 

if (edgeib .gt. O.dO) then 

c lm(f)=0 goes through the 00 to 01 edge (edge 3). 

I out =3 
else 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

lout=4 
end if 
else 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 

if (edgeib .It. O.dO) then 

c !m(f)=0 also runs through 00 to 10 and 00 to 01 edges, 

c Store crossing location and in/out information. 

knot12=knot12+1 
c loc12r(knot12)=lr 

c loc12i (knot 12)= l i 

end if 
end if 
c ******* 

go to 60 

c enter mesh square from top side or exit rectangle at bottom edge. 

56 l i = l i - 1 

if (li .ge. jbot) go to 57 
nre2=nre2+1 
kedge2(nre2)=lr+1 
go to 10 

57 fx01=fx00 
fxOI r=fx00r 
fxOI i=fx00i 
fxl 1=fxl0 
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f xl 1 r=f xIOr 
fxl 1 i=fx10i 
58 continue 

call f indfx( lr, l i , fxOO, fxOOr, fxOOi ) 
call f indfx( lr+1 , l i , fxlO, fxl Or, fxl 0 i ) 

Q* ★ ★ *★ ★ ★ 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei l = f xOO i *fx01 i 
edgei r=fx10i*fx1 1 i 
if (edgei l .gt. O.dO) then 
c lm(f)=0 goes through the 00 to 11 line, 

if (edgeir .gt. O.dO) then 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

l out =4 
else 

c lm(f)=0 goes through the 10 to 11 edge (edge 1). 

lout=1 
end i f 
else 

c lm(f)=0 goes through the 00 to 01 edge (edge 3) 

lout=3 

if (edgeir .It. O.dO) then 

c lm(f)=0 also runs through 00 to 10 and 10 to 11 edges, 

c Store crossing location and in/out information. 

knot23=knot23+1 
c loc23r(knot23)=l r 

c loc23i(knot23)=li 

end if 
end i f 
c 

£*★★★*★* 

60 continue 

nrsqu=nrsqu+1 

if(nrsqu .gt. maxnsq) go to 95 

c Test for there being at least one re(f)=0 line entering and 
c leaving the mesh square, 

c 

if ( (fx00r*fx10r .gt. O.dO) .and. (fx01r*fx11r .gt. O.dO) 

+ .and. (fx00r*fx01r .gt. O.dO)) go to (41,46,51,56) lout 
c 

c Computate the values of the modal function at the corners of a 
c a mesh square to determine its Taylor series to the 3rd order 
c for estimating its root locations, 
c 
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397 

398 

399 
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c f 00=one 

f 10=cdexp(fx10-fx00)-one 
f01=cdexp(fx01 -fx00)-one 
f 1 1=cdexp(fx1 1 - f xOO) -one 
c 

c * ****************************************** ****** ********************* 
c write (16,3001) ni,nf ,lr, li ,knot12,knot23, knot34,knot41 

c 3001 formate/ 1 ni , nf, Ir, ti and knot12, 23, 34 and 43 before ROOTS 
c + 2i6,2x,2i6,2x,4i6) 

c 

c *********** estimate locations of zeroes by radicals ***************** 
c 

call roots(f 10, f 01 , f 1 1 , sol ,nrsol ) 
c 

do 63 n=1,nrsol 

ureal = dreaKsol (n) ) 
uimag = dimag(sol(n)) 

if (ureal .It. O.dO .or. ureal .gt. I.OdO) go to 63 
if (uimag .It. O.dO .or. uimag .gt. I.OdO) go to 63 

62 thetaO )=( Ir+ureal )*tmesh 
theta(2)=( l i+uimag)*tmesh 
nf = nf+1 

zeros(nf )=dcmplx( theta( 1 ), theta (2)) 
nrzl=nrzl+1 - 

63 cont i nue 

c* ********************************************************** ******** 

c write (16,3002) ni,nf,nrsol 

c 3002 formate/ 1 out of ROOTS at 63, ni , nf and # of roots 1 , 3 i 4 ) 

c************************************************** ***************** 

c continue following the phase line 

go to (41,46,51,56) lout 
c ****** 
cc 

80 continue 

c 

return 

£***** 

95 continue 

write(16,9500) 

wri te( 16,4001 )lr, l i , ni ,nf ,tmesh 
wri te(* ,9500) 

4001 formate'go to 5 from 95 at Ir, li =' , i6, 1 , 1 , i6, ' ni, nf =',i6, 
♦ , ,',i6, , f mesh size =',d14.6) 
go to 5 

£***** 
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442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 

456 

457 

458 



97 continue 

wri te(16 # 9700) 

wri te( 16,4002) l r, l i # ni ,nf , tmesh 
write(* ,9700) 

4002 format ('go to 5 from 97 at Ir, li = 1 , i6, 1 , ' , i6, 1 ni, nf =',i6, 
+ , ,',i6 f ', mesh size =' ,d14. 6, /'zeroes found are kept.') 
c nf=ni 



c**** format statements 

9500 format(/5x, 1 too many squares on same phase line -- 1 , 

$ 'reduce tmesh and start over') 

9700 format(/5x, ' tmesh has been reduced but problems remain in', 
S 1 executing fzerox') 

c 

end 
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APPENDIX C: SUBROUTINE ROOTS 



This Appendix contains the listing of the subroutine ROOTS. This subroutine 
replaces the portion of the subroutine FZEROX where the coefficients of a quadratic 
equation are determined, and the subroutine QUAD for locating the zeroes of a 
quadratic polynomial. In the revised subroutine FZEROX, the roots of a cubic 
polynomial has to be found. This subroutine determines these zeroes by radicals. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 



subroutine roots (f 1 ,f 2, f3,sol ,nrsol ) 
£*★★****★★★★★**★*★*************★★★ ************************************* 
c This subroutine finds the roots of a third order polynomial by 
c radicals when the values of this polynomial at z=0, z=1 # z=i and 
c z=1+i are given as f 0=1 # fl+fO, f2+f0 and f3+f0 respectively, 
c Note that this algorithm takes cubic roots of two complex numbers 
c (hence the name 'solution by radicals') and use their linear 
c combinations as the roots of a third order polynomial. 
c ** *********************************************** ********************* 
implicit real*8 (a-h, o-z) 

complex* 16 fl , f2,f3, zero, one, ci ,ep14,em14,ep23,em23, 

+ fa, fb, fc,fd, fal , fa2, fa3, fa1s,p,q,delt, z,zm,u,v, sol 

parameter (xbi t52=52.d0*0.69314718055994531d0, thrd=1 .d0/3.d0, 

+ bi 1 50= 1 . d0/33554432. d0/33554432. d0,bit51=bi t50/2.d0, 

+ bi t52=bi t 51 /2.d0, tol =0.001 dO, 

+ zero=(0.d0,0.d0),one=(1 .d0 # 0.d0) ,ci=(0.d0, 1 .d0) , 

+ ep14=(0.5d0 # 0.5d0) ,em14=(0.5d0, -0.5d0) , 

+ ep23= ( - 0 . 5d0 , 0 . 86602540378443864675d0 ) , 

+ em23= ( - 0 . 5d0 , - 0 . 86602540378443864675d0 ) ) 

dimension sol(*) 
fa=one 

fb=(f2-ci*f 1+em14*f3) 
fc=((ep14+one)*f 1 - (em14+one)*f2+ci*f3) 
fd=(em14*(f2-f 1 )-ep14*f3) 
if (cdabs(fb) ,le. bit50) fb=zero 
if (cdabs(fc) .le. b i 1 5 1 > fc=zero 
if (cdabs(fd) .le. b i 1 5 2 ) fd=zero 
if (fd .ne. zero) then 
f a1=(- thrd)*fc/fd 
f a2=fb/fd 
f a3=fa/fd 
fa1s=f a1*fa1 
p=thrd*fa2-fa1s 

q=0.5d0*(f a3+f a1*fa2)-f a1*fa1s 
if (p .eq. zero) then 
if (q. eq. zero) then 
nrsol=1 
sold )=fa1 
return 
else 

nrsol=3 

u=((-2.d0)*q)**thrd 
sold )=u+fa1 
sol(2)=ep23*u+fa1 
sol(3)=em23*u+fa1 
return 
end if 
else 
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49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 



if (q. eq. zero) then 
nrsol=3 
sol(1)=fa1 
u=cdsqrt(( -3.d0)*p) 
sol(2)=fa1+u 
sol(3)=fa1-u 
return 
else 
v=p/q 
z=p*v*v 
absz=cdabs(z) 
if (absz .It. tol) then 
zm=-z 

fn=dint(1 .dO-xbi t52/dlog(absz)) 
lastn=idint(fn)-1 
dnn=fn-0.5d0 
dnd=fn+1 .OdO 
delt=one 

do 100 nt=1 , lastn 
dnn=dnn-1 .d0 
dnd-dnd- 1 .d0 

del t=(dnn/dnd)*delt*zm+one 
continue 

del t=(0.5d0*del t/q)**thrd 
u=p*delt 
v=-1 .d0/del t 
else 

de 1 1 =cdsqr t ( one+ z ) - one 
u=(q*delt )**thrd 

V s “ p/ u 

end if 
nrsol=3 
sol(1 )=u+v+fa1 
sol(2)=ep23*u+em23*v+fa1 
sol(3)=em23*u+ep23*v+fa1 
return 
end if 
end if 

else if (fc .ne. zero) then 
if (fb .eq. zero) then 
if (fa .eq. zero) then 
nrsol=1 
sold )=zero 
return 
else 

nrsol=2 

z=cdsqrt(-fa/fc) 
sol(1)=z 
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97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 
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121 
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127 

128 

129 

130 

131 

132 

133 

134 

135 
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138 



sol ( 2 )=- 2 
return 
end if 
else 

f al =0.5d0*fb/fc 

fa2=fa/fc 

Z=fa2/fa1/fa1 

absz=cdabs(z) 

if (absz .It. tol) then 

fn=dint(1 .dO-xbi t52/dlog(absz) ) 

lastn=idint(fn)-1 

dnn=fn-0.5d0 

dnd=fn+1 .OdO 

delt=one 

do 200 nt=1 , lastn 
dnn=dnn- 1 .d0 
dnd=dnd-1 .d0 

del t=(dnn/dnd)*del t*z+one 
continue 

delt=-0.5d0*del t/f al 
nrsol=2 

sold )=fa2*del t 
sol(2)=1 .d0/del t 
return 
else 

del t=cdsqrt(one-z) 
nrsol=2 

sol (1 )=-fa1*( one- del t) 
sol (2) = -fa1*(one+del t ) 
return 
end if 
end if 

else if (fb .ne. zero) then 
nrsol=1 
sol( 1 )=- fa/f b 
return 
else 

nrsol=1 
sold )=ep14 
return 
end if 
end 
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APPENDIX D: SUBROUTINE ABCOEF 



This Appendix contains the listing of the subroutine ABCOEF. The consistency 
self-checking procedure has been implemented to determine the correct method to 
evaluate the A - t and coefficients. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 



subroutine abcoef (zero,m) 



£***** 

c For each mode m, this suboutine calculates A-B coefficients in 

c all layers for combining two linearly independent solutions of 

c Stokes' equation to form the height gain function: 
c 

c height gain=exp(bcoefx( l ,m) )*(k1*exp(acoefx( l ,m) )+k2) 

c 

c where kl and k2 are two independent solutions to Stokes' 

c equation. In the top layer (i.e. nzlayr) the height gain is: 

c 

c height gain=exp(bcoefx( l ,m) )*h2 

c where h2 is a solution to the Stokes' equation associated 
c with outgoing energy flow. Here kl and k2 are proportional 

c to the kl and k2 used by Marcus and the h2 is proportional 

c to a modified Hankie function of order 1/3. 

c inputs... 

c zero-an eigenvalue in qll space 

c outputs... 

c acoefx-two dimensional array of complex exponents 

c coefficients used to combine two linearly 

c independent solutions of stokes' equation 

c bcoefx-two dimensional array of complex exponents 

c coefficients used for normalizing the height gains 

c note: acoefx and bcoefx are passed by the 

c common block /pap2/ 

c subroutines called... 

c xcdai 

c xcadd 

c common block areas... 

c coml 

c com2 

c papl 

c pap2 

£****★ 



implicit real*8(a-h,o-z) 

complex* 16 acoefx, bcoefx ,cqi j , h2xq1 ,dh2xq1 , h2xq2,dh2xq2 ( klxql , 
$ dklxql , k1xq2,dk1xq2 # k2xq1 ,dk2xq1 , k2xq2,dk2xq2,h2dk1x # 
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46 $ dh2k1x,h2dk2x,dh2k2x,numax,denax,numbx,denbx, intlx, int2x, 

47 $ hyx, dhyx,k1dhyx,dk1hyx,dk2hyx,k2dhyx, gamma, dgamdq, i , 

48 $ koa123,rtsumx,zero,q1 ,q2,sumx,surfno,dqi j,dqi jdz,sqng, 

49 $ dnunbx,dhux,dhlx,e13x,cneg,cldqzl ,cldqzm,cigama,koawav,tthd, 

50 + tacoef ,dacoef 

51 

52 parameter(downi=1 .d-3,downr=1 .d-3/0.4342944819032518d0, 

53 ♦ pi =3 . 1 4 1 592653589793238462643d0 , 

54 + i=(0.0d0,1.0d0),tthd=(2.d0/3.d0)*i, 

55 + cneg= ( 0 . OdO , 3 . 1 4 1 592653589793238462643d0 ) , el 3x=cneg/3 . dO ) 

56 c***** 

57 c mx l ayr=max i mum number of layers allowed 

58 c mxmode=maximum number of modes allowed 

59 

60 c 

61 c use include file for parameters of 

62 c use include file for parameters of 

63 c mxlayr max # layers 

64 c mxmode max # modes 

65 c 

66 Sinclude: 'mlaparm. inc 1 

***** Begin listing of: mlaparm. inc 

1 c 

2 c include file to define the 

3 c maximum # of layers (mxlayr) 

4 c maximum # of modes (mxmode) 

5 c 

6 parameter (mxlayr=35 ) 

7 parameter (mxmode=127) 

***** End listing of: mlaparm. inc 

67 c 

68 c 
£****★ 

70 c acoefx-two dimensional complex array used for combining two 

71 c independent solutions to stokes' equation 

72 c bcoefx-two dimensional complex array used for normalizing height 

73 c gain 

74 c cqij-two dimensional array containing coefficients for evaluating 

75 c qij in terms of qll 

76 c dqij -array containing coefficients for evaluating qij in terms of 

77 c qll 

78 c dqijdz-array containing derivatives of qi(z) in the different 

79 c layers 

80 c zi -array containing input hesights for the modified refractivity 

81 
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82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 



dimension acoefx(mxlayr ,mxmode) , 

$ bcoefx(mxlayr,mxmode), 

$ dqi j(mxlayr),cqi j(mxlayr f 2), dqi jdz(mxlayr) f zi (mxlayr+1 ) 



common /coml/f req,waveno,sqng 
common /com2/cqi j # dqi j ,dqi jdz,nzlayr 
common /papl/nrmode, koa123,surfno,zi 
common /pap2/acoefx f bcoefx 



£*★*** 

c check for single layer 

c 

c set a complex variable koawav=- i*koa123/(waveno*waveno) to 
c avoid repeating computations 

koawav=- i*koa123/(waveno*waveno) 

if(nzlayr .eq. 1)then 
q1=cqi j(1,1 )+zero*dqi j ( 1 ) 
call surf (ql f gamma, dgamdq) 

call xcdai ( -ql ,k2xq1 ,dk2xq1 , klxql ,dk1xq1 ,h2xq1 ,dh2xq1 ) 
dh2xq1=dh2xq1+e13x 

int1x=cdlog(koawav*dgamdq-q1/dqi jdz(1 ) )+2.0d0*h2xq1 
int2x=2.0d0*dh2xq1 -cdlog(-dqi j dz ( 1 ); 
call xcadd(sumx, intlx, int2x) 
rtsumx=0.5d0*sumx 
bcoefx(1 ,m)=-rtsumx 
return 
end if 

cldqzl=cdlog(-dqi jdz(1 )) 

c if l equals one then initialize cumulants and caculate a’s and 
c b's in bottom layer using ground boundary conditions. 

ql =cqi j ( 1 , 1 )+zero*dqi j ( 1 ) 

call xcdai (-ql # k2xq1 ,dk2xq1 , klxql # dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1 -e13x 

call surf (ql f gamma, dgamdq) 

c i gama=cd l og( i *gamma ) 

cal l xcadd(numax,cldqzl -cneg+dk2xq1 ,ci gama+cneg+k2xq1 ) 
cal l xcadd(denax,cigama+ klxql ,cldqzl+dk1xq1 ) 
acoefx(1 ,m)=numax-denax 
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127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 



call xcadd(denbx,k2xq1 ,acoefx( 1 ,m)+k1xq1 ) 
bcoefxO ,m)=-denbx 

calculate contributions to normalizing integrals. 

cal l xcadd(hyx,k2xq1 ,acoefx(1 ,m)+k1xq1 ) 
hyx=bcoef x( 1 ,m)+hyx 

cal l xcadd(dhyx,dk2xq1 ,acoefx(1 ,m)+dk1xq1 ) 
dhyx=bcoefx(1 ,m)+dhyx 

int1x=cdlog(koawav*dgamdq-q1/dqi jdz ( 1 ))+2.0d0*hyx 

int2x=2.0d0*dhyx*cldqzl 

call xcadd(sumx, intlx, int2x) 

do 9010 l=2,nzlayr-1 
lm1=l -1 

cldqzl=cdlog(-dqi jdz(l)) 
cldqzm=cdlog(dqi jdz(lml)) 
q1=cqi j( 1,1 )+zero*dqi j(l) 

cal l xcdai ( *q1 , k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1) 

dk2xq1 =dk2xq1 +cneg 

dk1xq1=dk1xq1 -e13x 

q2=cqi j( 1ml ,2)+zero*dqi j( 1ml ) 

cal l xcdai (-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

cal l xcadd(hyx,k2xq2, acoef x( 1ml ,m)+k1xq2) 

cal l xcadd(dhyx,dk2xq2,acoefx( 1ml ,m)+dk1xq2) 

k1dhyx=k1xq1+dhyx 

dk1hyx=dk1xq1+hyx 

dk2hyx=dk2xq1 +hyx 

k2dhyx=k2xq1+dhyx 

call xcadd(denax,cldqziTH-k1dhyx,cldqzl+dk1hyx) 

ca l l xcadd( numax, c Idqz l -cneg+dk2hyx, c ldqzm+cneg+k2dhyx) 

acoef x( l ,m)=numax-denax 

call xcadd(denbx,k2xq1 , acoef x( l ,m)+k1xq1 ) 

numbx=bcoef x( 1ml ,m)+hyx 

dnumbx=bcoefx( Iml ,m)+dhyx 

bcoefx( l ,m)=numbx-denbx 

calculate contribution to normalizing integrals. 

int1x=cdlog(-q1/dqi jdz( l )+q2/dqi jdz( Iml ))+2.0d0*numbx 
call xcadd(sumx,sunx, intlx) 
cal l xcadd(dhux,dk2xq1 ,acoefx( l ,m)+dk1xq1 ) 
dhux=bcoef x( l ,m)+dhux 
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172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 



i nt 1 x=2 . 0d0*dnumbx - c l dqzm 
int2x=2. OdO*dhux-c Idqzl 
call x c add (s unix ,sumx, intlx) 
call xcadd(sumx,sumx, int2x) 
9010 continue 



c if l equals nzlayer, calculate a's and b's using outgoing 
c wave in top layer. 

nzm1=nzlayr-1 

q1=cqi j(nzlayr, 1 )+zero*dqi j(nzlayr) 

call xcdai (-ql , k2xq1 ,dk2xq1 ,k1xq1 f dk1xq1 ,h2xq1 # dh2xq1 ) 

dh2xq1 =dh2xq1+e13x 

q2=cqi j(nzm1 ,2)+zero*dqi j(nzml) 

ca 1 1 xcda i ( - q2 , k2xq2 , dk2xq2 , k 1 xq2 , dk 1 xq2 , h2xq2 , dh2xq2 ) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

call xcadd(hyx # k2xq2,acoefx(nzm1 ,m)+k1xq2) 
numbx=bcoefx(nzlayr-1 ,m)+hyx 
bcoef x ( nz l ayr , m) =numbx - h2xq1 

c calculate contribution to cumulants. 

int1x=cdlog(-q1/dqi jdz(nzlayr )+q2/dqi jdz(nzm1 ))+ 

S 2.0d0*numbx 

call xcadd(sumx # sumx, intlx) 

cal l xcadd(dhyx,dk2xq2,acoefx(nzm1 ,m)+dk1xq2) 

dnumbx=bcoefx(nzm1 ,m)+dhyx 

int1x=2.0d0*dnunbx-cdlog(dqi jdz(nzm1 ) ) 

call xcadd(sumx,sumx, intlx) 

dhux=bcoefx(nz Iayr # m)+dh2xq1 

i nt2x=2 . 0d0^dhux-cd log( -dqi jdz( nz l ayr ) ) 

call xcadd(sumx # sumx # int2x) 

c renormalize b’s so that height gain integral equals unity. 

rtsumx= . 5d0*sumx 
do 9000 ll=1,nzlayr 

bcoef x( l l ,m)=bcoefx( l l ,m)-rtsumx 
9000 continue 
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218 

219 
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227 
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234 

235 

236 

237 

238 

239 

240 

241 

242 
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l=nzlayr 

lm1=l-1 

cldqzm=cdlog(dqi jdz( Iml )) 
cldqzl=cdlog(-dqi jdz(l)) 

c calculate q and associated quantities at bottom of layer l 
q1=cqi j ( l , 1 )+zero*dqi j ( l ) 

call xcdai ( -ql , k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 f h2xq1 ,dh2xq1 ) 
dh2xq1=dh2xq1+e13x 

q2=cqi J ( Iml ,2)+zero*dqi j ( Iml ) 

call Xcdai(-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

c Caculate acoefx( Iml # m) # bcoefx( Iml # m) 

c and cumulants using outgoing wave in nzlayr 

dh2k1x=dh2xqUk1xq2 
h2dk1 x=h2xq1 +dk1 xq2 
h2dk2x=h2xq1 +dk2xq2 
dh2k2x=dh2xq1+ k2xq2 

call xcadd(denax,cldqzl-cneg+dh2k1x,cldqznH-cneg+h2dk1x) 
call xcadd(numax f cldqzm+h2dk2x,cldqzl+dh2k2x) 

c If in the nzlayr-1 layer the magnitudes of A coefficients from 
c integration up and down differ by less than 0.02 dB and their 
c phases differ by less than O.OOlpi, the A and B coefficients 
c obtained from integration up will be accepted. 

tacoef =numax-denax 
dacoef=tacoef -acoefxC Iml f m) 
dif r=dabs(dreal(dacoef )) 
if (difr .It. downr) then 
d i f i =d i mag ( dacoef ) /pi 
dif i=dabs(dif i-dnint(dif i/2.d0)*2.d0) 
if (difi .It. downi) return 
end if 

acoefx( Iml ,m)=tacoef 

call xcadd(denbx,k2xq2,acoefx( Iml ,m)+k1xq2) 
bcoefx( Iml ,m)=h2xq1 -denbx 
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c calculate contributions to cunulants 

sunx=cdlog( -ql/dqi jdz( l )+q2/dqi jdz( Iml ))+2.0d0*h2xq1 

call xcadd(dh lx,dk2xq2,acoef x( Iml ,m)+dk1xq2) 

dhlx=bcoefx(lm1 ,m)+dhlx 

int1x=2.0d0*dh2xq1 -cldqzl 

call xcadd( intlx, sunx, intlx) 

int2x=2.0d0*dhlx-cldqzm 

call xcadd(sunx, intlx, int2x) 

do 9030 l=nzlayr- 1 ,2,-1 
lm1=l -1 

cldqzl=cdlog( -dqi jdz( l )) 
cldqzm=cdlog(dqi jdz( Iml )) 

c calculate q and associated quantities at bottom of layer l 
q1=cqi j( l , 1 )+zero*dqi j ( l ) 

call xcdai ( -ql ,k2xq1 f dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 , dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1 -e13x 

q2=cqi j( Iml ,2)+zero*dqi j( Iml ) 

call xcda i ( - q2 , k2xq2 , dk2xq2 , k 1 xq2 , ak 1 xq2 , h2xq2 , dh2xq2 ) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

dh2xq2=dh2xq2+e1 3x 

£***** 

c Calculate acoefx( Iml ,m) ,bcoefx( Iml # m) and cunulants 

c using continuity relations in terms of the linearly 

c independent functions kl and k2 

call xcadd(hyx,k2xq1 ,acoefx( l f m)+k1xq1 ) 

call xcadd(dhyx,dk2xq1 ,acoef x( l ,m)+dk1xq1 ) 

k 1 dhyx=k 1 xq2+dhyx 

dk1hyx=dk1xq2+hyx 

dk2hyx=dk2xq2+hyx 

k2dhyx=k2xq2+dhyx 

call xcadd(denax,cldqzl-cneg+k1dhyx f cldqzm+cneg+dklhyx) 
call xcadd(nunax f cldqzm+dk2hyx,cldqzl+k2dhyx) 
acoefx( Iml ,m)=nunax-denax 
call xcadd(denbx # k2xq2 f acoefxC Iml f m)+k1xq2) 
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numbx=bcoefx( l ,m)+hyx 
dnumbx=bcoefx( l ,m)+dhyx 
bcoefx( 1ml ,m)=numbx-denbx 

c calculate contributions to cunulants. 

int 1x=cdlog(-q1/dqi jdz( l )+q2/dqi jdz( 1ml ))+2.0d0*ntinbx 

call xcadd(sumx,sumx, intlx) 

cal l xcadd(dhlx,dk2xq2,acoefx( 1ml ,m)+dk1xq2) 

dhlx=bcoefx( 1ml ,m)+dhlx 

i nt 1 x=2 . 0d0*dnumbx- c l dqz l 

i nt 2x=2 . 0d0*dh l x- c l dqzm 

call xcadd(sumx,sumx, intlx) 

call xcadd(sumx,sumx, int2x) 

9030 continue 

c if l equal to one calculate ground 

c contribution to cumulants and renormalize bcoefx’s 

1 = 1 

q1=cqi j ( l , 1 )+zero*dqi j( l ) 

cal l xcdai ( -ql ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 
dk2xq1 =dk2xq1 +cneg 
dklxql =dk1xq1 -e13x 

call xcadd(hyx, k2xq1 ,acoefx( l ,m)+k1xq1 ) 
call xcadd(dhyx,dk2xq1 ,acoefx( l ,m)+dk1xq1 ) 
call surf (ql , gamma ,dgamdq) 
nixnbx=bcoefx( l # m)+hyx 
dnumbx=bcoefx( l ,m)+dhyx 

int1x=cdlog(koawav*dgamdq-q1/dqi jdz( l ))+2.0dO*numbx 
int2x=2.0d0*dnumbx-cdlog(-dqi j dz< 1 )) 
call xcadd(sumx # sumx, int lx) 
call xcadd(sumx,sumx, int2x) 

c renormalize b's so that height gain integrals equal unity. 

rtsumx= . 5d0*sumx 

do 9020 l 1=1 ,nzlayr-1 

bcoefx( l l ,m)=bcoefx( l l ,m)-rtsumx 

9020 continue 
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bcoefx(nzlayr,m) = - rtsimx 



352 


! 


353 




354 




355 


return 


356 


end 



67 



LIST OF REFERENCES 



1. D.E. Kerr, Propagation of Short Radio Waves, Peregrinus Ltd, London, United 
Kingdom, 1987. 

2. S.W. Marcus (1982), "A model to calculate EM fields in tropospheric duct 
environments at frequencies through SHF," Radio Science 17(5), 1108-1124 • 

3. V.I. Fock ( 1965), Electromagnetic Diffraction and Propagation Problems, 4 14 + ix 
pp., Pergamon Press, New York. 

4. L.W. Yeoh (1990), "An analysis of M-Layer: a multilayer tropospheric 
propagation program," Technical Report NPS-62-90-009, Naval Postgraduate 
School, Monterey, California 93943. 

5. D.G. Morfitt and C.H. Shellman, "MODESRCH, An improved computer 
program for obtaining ELF/VLF/LF mode constants in an earth-ionosphere 
waveguide," Interim Report 77T, Naval Ocean Systems Center, San Diego, CA 
92152, October 1976. 

6. Z. Schulten and D.G.M. Anderson, "An algorithm for the evaluation of the 
complex Airy functions", Journal of Computational Physics, Vol. 31, No. 60-75, 
1979. 



68 



INITIAL DISTRIBUTION LIST 



1. Defense Technical Information Center 1 

Cameron Station 

Alexandria, VA 22304-6145 

2. Library, Code 52 1 

Naval Postgraduate School 

Monterey, CA 93943-5002 

3. Department Chairman, Code EC 1 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey, CA 93943-5100 

4. Professor Hung-Mou Lee, Code EC/Lh 1 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey, CA 93943-5100 

5. Professor Lawrence J. Ziomek, Code EC/Zm 1 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey, CA 93943-5100 

6. Han, Yin Yuan 2 

5F1, No. 10, Lane 165, HSIN-LONG RD. SEC. 4, 

Taipei Taiwan, R.O.C. 

7. Naval Academy Library 1 

Kao-Hsiung Tso-Ying P.O. 90175 

Taiwan, R.O.C. 

8. Ting, Chueh 1 

Department of Electrical Engineering 

Chung Cheng Institute of Technology 
Tao-Yuan, Tai-Hsi 
Taiwan, R.O.C. 



69 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOl 
MONTEREY CA 93943-5101 



GAYLORD S 



